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Abstract. With the example of the spherically symmetric scalar wave equation on 
Minkowski space-time we demonstrate that a fully pseudospectral scheme (i.e. spectral 
with respect to both spatial and time directions) can be applied for solving hyperbolic 
equations. The calculations are carried out within the framework of conformally com- 
pactified space-times. In our formulation, the equation becomes singular at null infinity 
and yields regular boundary conditions there. In this manner it becomes possible to 
avoid "artificial" conditions at some numerical outer boundary at a finite distance. We 
obtain highly accurate numerical solutions possessing exponential spectral convergence, 
a feature known from solving elliptic PDEs with spectral methods. Our investigations 
are meant as a first step towards the goal of treating time evolution problems in General 
Relativity with spectral methods in space and time. 
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1. Introduction 

The simulation of dynamical processes within the theory of General Relativity plays 
an important role for understanding astrophysical phenomena, for studying the sta- 
bility of equilibrium configurations by introducing small perturbations and evolving 
in time and for predicting the properties of the emitted gravitational waves. An ideal 
way of analyzing such processes carefully would be the construction of explicit solu- 
tions to Einstein's field equations. However, due to the mathematical complexity of 
these equations, even stationary configurations can be described in terms of analytic 
expressions in only a few exceptional cases. Therefore, the only chance for tackling 
time dependent processes is the application of numerical methods. Nevertheless, 
in order to come as close as possible to an explicit solution, it is desirable to find 
mathematical descriptions and numerical procedures that permit the computation 
of a very accurate approximation to the solution in question. 

As we demonstrate in this paper, a promising approach towards this goal is 
the combination of highly accurate (pseudo-) spectral methods with the fruitful 
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concept of conformally compactified space-timecJ. By means of spectral methods, 
general relativistic equilibrium configurations have been obtained with almost ma- 
chine accuracy [3] . Dynamical relativistic problems have also been studied utilizing 
spectral methods (with respect to the spatial directions, combined with finite dif- 
ference methods in time direction), see e.g. [4]. (For a comprehensive overview of 
the applications of spectral methods in general relativity see [10].) However, there 
is only little experience regarding spectral expansions with respect to space and 
timet^ px3j . As a first step towards the goal of treating time evolution problems in 
General Relativity with a fully pseudospectral scheme, we study model equations, 
in particular linear and non-linear wave equations. 

The concept of conformal infinity^ i.e. conformally compactified space-times, was 
introduced by Penrose [14j in 1964 (For an overview on this topic see [7]; numerical 
studies can be found in |8|ll|12|r7] ). Within this scheme we are able to carry out 
the numerical calculations up to infinity. As a consequence, we have no need of 
numerical outer boundaries at a finite physical distance. At such finite boundaries 
one usually imposes particular conditions in order to complete the mathematical 
problem. These conditions have to be compatible with the differential equation 
to be solved and should lead to a well-posed problem. (Pseudospectral methods 
are particularly sensitive with respect to this issue. In the worst case, the numerical 
simulation breaks down after some time when the errors arising due to incompatible 
boundary conditions accumulate.) However, in any case the physical meaning of such 
boundary data and their influence to the solution is unclear. That is why we prefer 
the conformal compactification. 

An additional important feature of the conformal approach is the precise deter- 
mination of gravitational wave signals at null infinity .y , which is located at finite 
coordinate distance. Hence one avoids approximative wave extraction techniques at 
finite physical distance. 

Within the conformal concept, hyperbolic differential equations become singular 
at null infinity. Although it is possible to carry out an appropriate regularization in 
some special cases (for Einstein's field equations, a suitable reformulation in terms 
of Friedrich's regular conformal field equations [9; can be used), we demonstrate 
here that this degeneration of the equations permits a careful numerical treatment. 
In combination with a fully pseudospectral scheme, we obtain highly accurate nu- 
merical solutions (up to 12 or 13 correct digits for a double precision code). 

The paper is organized as follows. In Sec. [21 we recall the conformal compactifi- 
cation of Minkowski space-time. Moreover we consider the scalar wave equation on 
this background and derive the boundary conditions to be imposed at the singular 
points of this equation. The numerical method for solving the singular equation 

''In principle, it is also possible to apply the pseudospectral algorithm described in this paper to 
arbitrary compact domains if the boundary conditions at are replaced by other conditions at 
the boundaries of these domains. 

''Only in the context of finite and spectral element methods the simultaneous space- time treatment 
of hyperbolic equations is already more common, see e.g. |15| and references therein. 
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is explained in detail in Sec. [3l In the subsequent section, we give a number of 
numerical examples in order to test the method and to study the accuracy of the 
numerical solutions. In Sec. El we show that the application of the spectral scheme 
is not restricted to the homogeneous linear wave equation. To this end we study two 
additional example equations: an inhomogeneous wave equation and a non-linear 
wave equation. In Sec.[6l we consider a regularized version of the wave equation and 
demonstrate that our numerical method is also applicable to this equation. Finally, 
in Sec. [7] we discuss our results. 



2. Wave Equation and Compactification 
2.1. The wave equation 

The model equation to be studied throughout most of this paper is the spherically 
symmetric wave equation 

□/ - f,rr + -f.r - f,tt =0, / = /(r, t) (2.1) 

r 

on Minkowski space-time with the line element 

ds^ = dr^ + r^{di}^ + sin2??d(^2) - dt'^ (2.2) 
using spherical coordinates {r,{},ip,t). The general solution to (|2.ip . 

fir.t)^^±±^l±^^, (2.3) 
r 

with a,b G C^(R), can be used to investigate the accuracy of the numerical solution. 
Here we restrict ourselves to solutions which are regular at r = 0. As a consequence, 
we obtain the condition b{x) = — a(— x) for all a; G R, i.e. we consider solutions of 
the form 

(2.4) 

For the conformal compactification of Minkowski space-time, one often uses the 
following coordinate transformation (see [7]): 

r = i [tan(i? + T)+ tan(i?, - T)] , t=^ [tan(i? + T) - tan(i? - T)] , (2.5) 

through which new coordinates R and T are introduced, with i? G [0,^], T G 
[R — ^, ^ ~ R]. The line element takes the form 

ds^ = [4(di?2 - dT^) + sin2(2i?)(di?2 + sin^d dip"^)] (2.6) 

with 

Q -.^ V2cos{R + T), P -.^ V2cos{R~T). (2.7) 

Obviously, the metric components are singular at null infinity .y , i.e. at R±T — tt/2 
where P = or Q = 0. However, the rescaled line element ds^ — Q^P^ds^ is regular 
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Fig. 1. The standard conformal diagram of Minkowski space-time. The Hne R = corresponds 
to the center of symmetry r = 0. Furthermore, past and future null infinity past and future 

time-like infinity , and space-like infinity are marked accordingly. 



at . The domain on which R and T are defined is illustrated by the well-known 
standard conformal diagram, see Fig. [TJ 

For the wave equation in terms of the coordinates R and T one obtains 

QPilRR - Itt) + 2 0. (2.8) 

The general solution (|2.4p can be written as 

f(R T) = MT + R)'A{T-R) 

■'^ ' ' tan(i? -f T) + tan(i? - T) ^ ' ' 

with A{x) :— a(tanx) for all x E [—tt/2, 7r/2]. Due to the appearance of the functions 
P, Q and sin(2i?), the wave equation (|2.8p is singular at infinity (as a consequence 
of the compactification) and at the center of symmetry (as an effect of the spherical 
coordinates), i.e. at the following points: 

• r and J^f~: i?-r= f P = 

• i+ and y+: R + T = ^ => Q = 

• i^: R±T = ^ P = Q = 

• i? = : sin(2i?) = 0. 

In the numerical scheme we need to treat these points particularly carefully, in 
order to derive well-defined boundary conditions, as will be discussed in the next 
subsection. 

We study two different types of initial boundary value problems (IBVP) for the 
wave equation (|2.8p . see Fig. O In a hyperboloidal IBVP we prescribe initial data 
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(the values of / and the time derivative /_t) on the hyperboloidal sliced T — Tmin 
and evolve these data up to another hyperboloidal slice T = T,„ax- For the second 
type we consider a "standard" Cauchy problem, in which initial data are given on 
the particular Cauchy surface T — t — 0. From these data, the entire future of / up 
to and ^+ is determined. As a consequence of the time symmetry of the wave 
equation, the past of / can be calculated in the same manner. We thus obtain the 
function / everywhere. 

For relativistic time evolution problems, initial data on Cauchy surfaces have 
to be constructed very carefully in order to avoid (logarithmic) singularities at i° 
which cannot be removed easily by means of a coordinate transformation. Solutions 
to this problem have been discussed by Corvino [5J, utilizing a "gluing" technique 
in order to obtain initial data which are exact Schwarzschild data near spatial 
infinity, and more generally by Dain and Friedrich [6]. Note that, on the other 
hand, there are specific conditions which guarantee the regularity at of initial 
data on hyperboloidal slices |1|2) . For this reason, we concentrate primarily on the 
hyperboloidal IB VP problem here. 




Fig. 2. The numerical domains for the hyperboloidal initial boundary value problem (left) and the 
standard Cauchy problem (right). 



2.2. Boundary conditions 

The appropriate boundary condition^ at the singular boundaries ^ and R = 
follow from an analysis of the wave equation (|2.8p at these "critical" points. In this 

"^A slice is called hyperboloidal if it is space-like everywhere and extends up to .y. 
"^Note that we denote any equations to be imposed at inner and outer borders of the numeri- 
cal domain as "boundary conditions", even if they may be of quite different nature (regularity 
conditions at coordinate singularities, conditions at , second order conditions, . . . ). 
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investigation we assume regularity of the solution, i.e. bounded first and second 
order derivatives. 

In order to derive a condition at i? = 0, we first multiply (|2.8p with sin(2i?) 
and then perform the limit R 0. Since P = Q — cosT as i? ^ 0, the following 
Neumann condition 

f.R = (2.10) 

arises. Being a necessary requirement in the context of spherical symmetry, this 
condition also appears as a consequence of the wave equation. 

At J^'^ the relation Q = holds, and hence the wave equation implies the 
condition 

Im - f.T = 0. (2.11) 

Similarly, at one finds that f^R + f^x = 0. Thus the tangential derivative of / 
along vanishes, i.e. f,fj,t^ = 0, where t^^ = (1, i'^, t'^, =Fl) is a tangential vector on 
J^. As a consequence, one may choose either the boundary condition (|2.1ip at 
or the Dirichlet condition 

/ = constant. (2-12) 

In the latter case, the constant needs to be read off from the initial data which also 
extend up to y. 



3. Numerical Method 

In this section we describe the numerical procedure for solving the wave equation 
([2?8ll with the boundary conditions (|2.10p and (|2.11|) or (|2.12|1 for the two types of 
mathematical problems — the hyperboloidal IBVP and the Cauchy problem (see 
Fig. ED. 

The numerical method consists of the following ingredients: 

(1) Mapping of the physical domain onto a unit square (or onto several unit squares) 
by introducing appropriate coordinate transformations. The coordinates being 
defined on the unit square(s) are referred to as spectral coordinates. 

(2) Expressing the "wave function" / in terms of another unknown function such 
that the initial conditions are satisfied automatically. 

(3) Expansion of this unknown function in terms of a truncated series of Chebyshev 
polynomials with respect to the spectral coordinates. A particular finite resolu- 
tion for the numerical approximation is chosen and appropriate grid points in 
the spectral coordinates on the unit square(s) are identified. 

(4) Formulation of an algebraic system of equations for the values of the unknown 
function at the coordinate grid points. This system results from the evaluation 
of the wave equation and the boundary conditions within the spectral approx- 
imation scheme being chosen. 
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(5) Calculation of the solution of this system by means of the Newton- Raphson 
method. 

These points are described in some detail in the following subsections. 



3.1. Spectral coordinates 

Smooth functions can be expressed in terms of spectral Chebyshev expansions pro- 
vided these functions are defined on an interval or, for more-dimensional functions, 
on a cross product of intervals. For this reason we introduce an appropriate coordi- 
nate mapping through which the physical domain is obtained as the image of the 
unit square (in the l-|-l-dimcnsional problem considered here). 

For the hyperboloidal IBVP we use the coordinate transformation 



R = 



^miii)'^ 



CT, T = T„; 



(3.1) 



with (T,T E [0,1]. An illustration of these coordinates is given in Fig. [5^, where 
particular a- and r-coordinate lines are shown. The physical boundaries T = Tmin, 
T = Tmax, R = and ^+ are mapped to the edges t — 0,t=1, a = and (7=1 
of the square, respectively. 

Note that at the points z° and the wave function is not analytic with respect 
to the coordinates R and T. Therefore, in the case of the standard Cauchy problem, 
it is necessary to introduce a singular coordinate mapping through which the wave 
function becomes analytic with respect to the spectral coordinates. To this end we 
divide up the physical domain into two triangular subdomains and map each of 
these subdomains onto a unit square, see Fig. [8)3. We use the particular coordinates 



R^ - [2a+{l-a)T] 



T= -(1-<j)t 



in domain 1 and 



R 



fT)r, T 



(1 - a)T] 



(3.2) 



(3.3) 



0, i", i+, 



in domain 2. The edges of the two squares correspond to i? = 0, T 
sections of J^+, and, additionally, the common boundary B = {(i?, T) : R = T;Q < 
R < between the two domains. Boundary conditions at B can be derived 

from the requirement that the wave function / be analytic in a neighbourhood of 
the initial slice T = 0. As a consequence of the wave equation, / is then analytic 
everywhere in the domain of dependence. In particular, for the standard Cauchy 
problem it follows that / possesses a continuous and differentiable transition at B. 

The particular choice of the spectral coordinates is not unique. Also the bound- 
ary B (in the Cauchy problem) could have been chosen to be at a different location 
(e.g. at T = aR with a — constant 7^ 1). Although our choice (T = R) represents 
a characteristic curve of the wave equation, in general the curve does not need to 
possess specific features with respect to the underlying equation. 
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(a) Spectral coordinates — Hyperboloidal initial value problem 



R = 




R = 



R = 




(b) Spectral coordinates — Cauchy problem 
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Fig. 3. Illustration of the spectral coordinates a and t for the two types of IBVPs. 
3.2. Initial conditions 

A prescribed set of initial conditions for the wave equation, to be imposed at T = 
Train (i-c at T = 0), Can be satisfied automatically through a specific ansatz for the 
wave function /. For given functions G{R) and H{R) with 



f{R, T = T^in) = G{R), It{R, T = T^in) = H{R), 



(3.4) 



we replace /, expressed as a function of the spectral coordinates a and r, by a 
function /2 via 



f{a,T)=fo{a)+h{a)T + h{a,T)T\ 



(3.5) 
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The functions fo{(j) and /i(cr) are given in terms of the initial data by 
/o(a) = /(a,T = 0) = G[i?(a,r = 0)] 

and 



/i(^) 



d£dR d£dT 



dGdR 



OT 



(3.6) 



(3.7) 



In this manner, for any regular choice of the new unknown quantity /2(ct, r), the 
function /(cr, r) as given in (j3.5p satisfies the correct initial conditions at r = 0. 



3.3. Chebyshev approximation 

In a spectral method the unknown functions are approximated in terms of appropri- 
ate basis functions. Here we choose Chebyshev polynomials defined on the interval 
[0, 1]. In particular, we approximate /2 by 

h{a,T) « ^^c,,r,(2a- l)r,(2r- 1), (3.8) 

i=0 j=0 

where T„ denotes the nth Chebyshev polynomial, Cij the Chebyshev coefficients^ 
and 

n^=N^ + l, nr = Nr + l (3.9) 

are the prescribed resolution orders with respect to the spectral directions. 

In accordance with this choice we introduce the following spectral collocation 
points Pij = {ai, Tj) at which the the wave equation and the boundary conditions are 
evaluated in order to build up an algebraic system of equations (see next subsection): 

i^O,l,...,N, (3.10) 
J = 0,l,...,iV,. (3.11) 

Note that we choose to use the extrema of the Chebyshev polynomials (Gauss- 
Lobatto collocation points) so as to have gridpoints lying on the boundaries. 



CT, — sm 

2 / j 

sm 

2Nr 



3.4. Algebraic system of equations 

For a given spectral approximation order, it is straightforward to compute from the 
values X^j /2(cri,rj) 

(1) the Chebyshev coefficients of /2 

(2) the Chebyshev coefficients of the first and second derivatives of f-z with respect 
to the spectral coordinates a and r 

(3) the values of these derivatives at the spectral collocation points Pij — {ui ,Tj). 



10 Jorg Hennig and Marcus Ansorg 



For any values Xij it thus becomes possible to "evaluate" the wave equation at the 
spectral gridpoints by inserting the function and derivative values. The mathemat- 
ical task to be solved can therefore be formulated as follows: Calculate the x n,- 
unknown values Xij as the solution of the algebraic system of the x rir equations 
Fij{Xki) — 0, where Fij is the left hand side of the wave equation or boundary 
condition evaluated from Xki at the grid-point {i,j). 

In the case of the hyperboloidal IBVP the following equations arise: 

at i? - : = 0, at y+ : f^R - /,t = 0, 

otherwise: the wave equation, (3-12) 

and in the Cauchy problem 

at R^O: f^R = 0, at ./+, io and i+ : f = 0, 
at B: f and ^ continuous, otherwise: the wave equation, (3.13) 

where df/dn denotes the normal derivative with respect to the boundary B (here, 
df /dn cx f,R — f^r)- The condition / = at infinity is true for all initial data 
vanishing at infinity. (/ = at infinity is automatically guaranteed for all bounded 
functions A{x) in the general solution (|2.9p of the wave equation.) 

It turns out that there are two exceptional points at which the conditions (|3.12|) 
or (|3.13p are already satisfied as a consequence of the "ansatz" p.Sp : 

Point 1: cr = r = 0, Point 2: cr = 1, r 0. 

These points correspond to the intersections of the initial slice T = Tmin with R — 
and with ^+ (in the hyperboloidal IBVP) or with i'^ (in domain 1 in the Cauchy 
problem). For example, at cr = r = 0, the condition f r = reduces to gR — 
which is already satisfied for the spherically symmetric initial data. 

Therefore, additional conditions need to be imposed at these two points, in order 
to complete the algebraic system of equations. As with the various other boundary 
conditions, these conditions also follow from the wave equation ()2.8p . 

In the limit R ^ one obtains 

3/,iii? -/,TT- 4 tanr/.T = 0. (3.14) 

At J^"*" the following condition can be derived 

fMR + f.TT - 21r,t = 0. (3.15) 

Finally, in the case of the Cauchy problem, at cr = 1 and r = we impose the 
condition 

/2 = 0. (3.16) 

This is a consequence of p.5p and the fact that / = at J^. Since /o = /i = 
for (T = 1, it follows that /2((t = 1, r 7^ 0) = 0. Hence, for a continuous function, 
/2(cr = l,r = 0) = holds. 
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3.5. Newton Raphson method 

Throughout most of this paper we concentrate on the hnear wave equation, for 
which the system Fij{Xki) = is hnear. A corresponding solution method as e.g. 
the Gauss- Jordan ehmination or the LU decomposition would provide us with the 
solution to our IB VP. However, in order to tackle non-linear equations it is necessary 
to use a more general scheme. For this reason we choose the Newton Raphson 
method. 

We define the x n^-dimensional vectors 

X := {Xqq.Xqi, . . . ,XoN^,Xio, . . . ,Xii\i^, . . . ,Xn^o, . . . ,Xiy^N^), (3.17) 
F := {FoQ,Foi . . . ,FoN^,Fio, . . . ,Fim^, ■ ■ ■ ,Fn^o, ■ ■ ■ ,Fn„n^), (3.18) 

containing all components of Xij and Fij. Hence, the system of equations can be 
written as 

F(X) = 0. (3.19) 

Within the Newton Raphson method, the system (I3.19|) is solved iteratively using 
an "initial guess" X° and computing the successive vectors 

X"+i = X" - [F'(X")]"iF(X"). (3.20) 

Here we approximate the Jacobi matrix F' = (dFi/dXj) through a second order 
finite difference representation and obtain the inverse via an LU decompositiorQ. 



4. Test of the Method with ExpUcit Examples 
4.1. Numerical accuracy 

We study several explicit examples in order to investigate the effectiveness of our 
numerical method. From (|2.9p we can construct a particular solution to the wave 
equation by choosing the appropriate function A{x). For a given time interval 
[Tmin, Tmax], whcrc < T^nin < Tmax < f (equal sigus in the case of the stan- 
dard Cauchy problem), we read off / and /t at T = Tmin as initial data for the 
numerical method and solve for / in the entire domain expanding up to the maximal 
time T = T„,ax. 

As a measure of the overall numerical accuracy, we calculate the numerical 
residual 

Rnumina^rir) ■■= max |/num(o-,T) - /an(cr, t)| , (4.1) 
(a,r)e[0,l]" 

where /num and /an denote the numerical and analytical solution, respectively. 

In order to investigate to what extent the numerical method is effective and 
leads to accurate solutions, we compare -Rnum to the analytical error caused by an 



'^Note that for our l-|-l-dimensional IBVP a direct matrix inversion is computationally feasible. 
However, for higher dimensional problems such methods are too computationally expensive and 
one requires the use of iterative matrix inversion methods. 
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exact Chebyshev representation of the given resolution order. This error is given by 
the foUowing analytical residual 

-Ran(^^o■,f^T) := max l/approxlo", r) - /an(a-, t)| (4.2) 
(cr,r)e[0,l]^ 

and describes the maximal difference between the analytical solution /an and its 
Chebyshev approximation 

n^ — l — 1 

/approx = J2 %-^.(2fT - 1)Tj{2t - 1) (4.3) 

of the order (ng-, n^- j^. 

Since the error (|4.2p of the Chebyshev approximation (|4.3p is close to the smallest 
possible polynomial approximation error (which one would obtain for the optimal 
approximation polynom of the particular function to be approximated), i?an pro- 
vides a good measure for the effectiveness and potential accuracy of our numerical 
method. In the subsequent subsection, we display plots of both residuals i?num and 
iian for several example solutions. 

A well-known property of Chebyshev approximations for smooth functions is an 
exponential decay of i?an with increasing resolution (up to a final saturation level of 
order 10""'^^ or 10""'^^ due to the finite machine accuracy of 16 numerical digits being 
used). As will be demonstrated below, the numerical residuals i?num also possess an 
exponential fall-off and reach a final saturation level near the final level of i?an- 

Another important issue is the dependence of the numerical residual on the 
size of the time interval [Tmin, Tmax]- It turns out that without substantial loss of 
accuracy we may choose a time step that is comparable to or larger than a spatial 
scale implied by the initial data. 



4.2. Numerical examples 

As a first example we consider A{x) — x in Eq. (12. 9|) . i.e. the analytical solution 

27? 

" tan(i? + T) + tan(i? - T) ' ^^'^^ 

(A plot of this solution can be found in Fig. [H]in | Appendix A[ ) The numerical and 
analytical residuals [as defined in (|4.ip . (|4.2p ] for a small and a large time interval 
in the hyperboloidal IBVP are shown in Fig. 01 We have chosen the resolutions Ua- 
and n-r to be connected by ria- = 2nT-. 

As expected, i?an exhibits geometric convergence (the graph is roughly a straight 
line in the logarithmic plot). For = 18 a saturation level is reached at a value 
below 10~^*. The numerical residual i?nmn shows the same qualitative behaviour: It 
decreases almost linearly in the logarithmic plot until it reaches a saturation level at 



For the calculation of the coefficients Cij we have used the exact values of the function /an at the 
gridpoints (ctj , Tj ) . 
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Fig. 4. Plot of the numerical and analytical residuals for the example solution with A{x) = x. 
Parameters: = 2n-r, Tmin = 0.3. Left panel: Tmax = 0.6. Right panel: Tmax = 1.2. (Note that 
with the latter value of Tmax a fairly large time interval is realized since future time-like infinity 
is characterized byT = 7r/2 = 1.57...) 



a value of about 10""'^ for tIt = 17 in the case of the small time interval and also at 
almost 10^^^ for rir — 20 in the case of the large time interval. Thus approximately 
the same accuracy of about 12 correct digits can be obtained for both the small 
and the large time interval. 

We study the solution in question with A{x) = x also as a standard Cauchy 
problem. The numerical result is presented in Fig. [51 The final saturation level 
of about 10^^^ (which is reached for 20 x 20 collocation points in each numerical 
domain) is even below the final level in the hyperboloidal problem. The reason for 
this might be the particular treatment of the points i° and i^, see Sec. 13.11 
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Fig. 5. Numerical accuracy in the standard Cauchy problem with A(x) = x. The numerical reso- 
lutions are chosen to be related by = «t in both subdomains. 



A number of further numerical examples for the hyperboloidal IBVP are sum- 
marized in Tab. [TJ (For plots of these solutions see Fig. [T^] in [Appendix A ) The 
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corresponding plots of the numerical accuracies are shown in Fig. [6l In these ex- 
amples, a saturation level between 10""'^° and 10"^'^ is reached for sufhciently large 
spectral resolutions. An exception is example F, which describes an initial pulse 
on a compact support. Even with a resolution of rio- = 160, rir = 40 (where the 
saturation level is not yet reached) the numerical solution is correct only in the first 
8 digits. This follows from the fact that a very accurate spectral approximation of a 
non-analytic C°°-function requires high resolution. As a consequence, the numerical 
techniques become expensive. 



Table 1. List of the further numerical examples for the hyperboloidal IBVP, see Figs.[6land ll2l 





A{x) 


Remark 


T ■ 

^ mm 


T 

max 




(A) 


sin(6a;) with 6 = 10 


arbitrary number of 
minima and maxima 
depending on b 


0.3 


0.6 


2nr 


(B) 


I sm{bx) with 6 = 20 


cf. (A) 


0.3 


0.6 


2nr 


(C) 




incoming "hill" 


0.3 


1.2 


2nr 


(D) 


g-64(2;-0.7)'^ 


incoming Gauss-like 
pulse, reflected at 
i? = 0, outgoing with 
inversed amplitude 


0.3 


1.0 


2nr 


(E) 


g-64(2:-l.l)^ _ g-64(2;+0.1)^ 


two Gauss-like pulses 
crossing each other 


0.3 


0.9 


2nr 


(F) 


1 

— 30e (=:+0.8){0.2-x) 

(if -0.8 < a; < 0.2, 
otherwise A{x) — 0) 


outgoing pulse with 
compact support 
{C°° , not analytic) 


0.3 


0.6 





5. Other Equations 

The numerical examples presented in the previous section provide evidence for the 
fact that highly accurate solutions to the wave equation can be obtained using 
a fully pseudospectral scheme. In this section we demonstrate that the method 
also works in the case of other differential equations. We discuss two particular 
modifications of the wave equation: an inhomogcneous wave equation and a non- 
linear wave equation. 
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Fig. 6. The residuals of the examples for the hyperboloidal IBVP in Tab.[T] see Fig. 1121 



5.1. Inhomogeneous wave equation 

An interesting feature of time evolution problems in General Relativity is the possi- 
ble formation of singularities from completely regular initial data, e.g. the scenario 
of the formation of a black hole from a collapsing star. For this reason we investigate 
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the applicability of our numerical method to an example solution which develops 
a pole like singularity. In particular we study how closely to this critical point the 
numerical domain can be located. Since solutions to the homogeneous linear wave 
equation do not develop singularities from regular initial data, we consider a par- 
ticular inhomogeneous linear wave equation, 

QP[f,BR.- f,TT) + 2 sm{2R) =nR,T), (5.1) 

where the right hand side I{R, T) is taken such that / develops a pole after a finite 
time. This property is guaranteed by choosing the solution 

/(^' P) = (i?2 _ ^2)2 „ y2)2 ' (5-2) 

which is singular for R ~ Rq, T = Tq. From this expression we can calculate 
the corresponding inhomogeneity / as well as the initial data. Equation (jS.ip can 
be solved in the same manner as the homogeneous wave equation. Note that the 
boundary conditions at and at the "exeptional points" , cf. Sec. 13. 4i take on a 
different form. At instead of (|2.1ip we impose the condition 

4sin(2i?)(/,fl,-/,T) =/. (5.3) 
For CT = r = we rewrite p.l4p accordingly: 

Sf.RR - f.TT - 4 tan Tf,T ^ „ (5.4) 

2 cos^ I 

Finally, at cr = 1, t = 0, we replace (|3.15p by the condition 

2 sin(2i?)(/,^;^^ + f^TT - 2f,RT) = I.R- (5.5) 

The residuals for the numerical solution of the hyperboloidal IBVP with the 
exact solution (|5.2|) are shown in Fig. [71 In [7K-c, the solution is calculated in a 
hyperboloidal domain with a fixed value of Tmin = 0.3 and increasing maximal time 
values Tiiiax approaching the singularity. One observes that more and more resolu- 
tion with respect to the time direction is required in order to reach the saturation 
level. Furthermore, the maximal accuracy obtained decreases slightly. However, even 
in a small vicinity of the pole, the solutions are very accurate for a sufficiently large 
resolution. 

In Fig. [7li, the inhomogeneous wave equation is solved within a narrow hyper- 
boloidal strip with Tmin = 1.1 and Tmax — 1-2. Although the singularity (located 
at i?o = 0.1, To = 1.3) is close to this domain, the solution is highly accurate 
(11 correct digits) and the numerical saturation is reached for a moderate value of 
rir = 20. 

The above examples demonstrate the applicability of our fully spectral method 
to solutions which encounter a singular behaviour. The numerical solutions being 
obtained retain high accuracy, provided that an appropriate resolution is chosen. 
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Fig. 7. The residuals Rnum and i?,an for solutions of the inhomogoneous wave equation, being 
calculated on different time intervals [Tj„in, Tmax]- The pole is located at Rq = 0.1, Tq = 1.3. A 
plot of the solution is shown in Fig. llOl in |Appendix A| 



5.2. Non-linear wave equation 

In the preceding sections we have studied linear differential equations. However, 
in view of future appUcations in General Relativity, it is interesting to apply our 
numerical method to non-linear equations. (In an appropriate formulation, the Ein- 
stein equations reduce to a set of non-linear wave equations.) 
To this end, we consider the example equation 



□/ ^ I, 



~f,r ^ f,tt — 

r 



1 



with A = constant. The general solution (regular at r = 0) is given b}Q 



(5.6) 



/(r-, t) = In [1 + A (a(i + r) - a(i - r))] 



(5.7) 



^The solution can be obtained by writing / as /(r, £) = <I>(r, t)/r and introducing null coordinates 
u = r-\-t^v = r — t. Then, Eq. 115.61 1 reduces to 3>,uu = which can be solved easily. 
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We introduce again the coordinates R and T [see (|2.5p ] and obtain the non-hnear 
wave equation 



QP 



QP{I, 



RR 



M + —7:7^ [iQ' + P')f,R + iQ' - P')f,T] 



sin(2i?) 



sin(2i?) " 



with Q and P as defined in (|2.7p . The general solution reads 

..^ ^^ ^ 2 \n[l + \{A{T + R)-A{T-K))] 
' > A tan(i? + T) + tan(i? - T) 



(5.8) 



(5.9) 



As with the linear wave equation, the boundary conditions can be obtained by 
analyzing the equation at the singular points. In this manner we find the conditions 



OV,fl = Af at i? = 0, /,«-/,T = at J+ 
and, additionally, at the "exceptional points" (see Sec. 13. 4p 

4A 



3/,_Rfl - /,TT - 4tanT/,T 



cos^T 



ff.R at = 0, r = 0, 



and 



f,RR + f,TT - 2/,flT = at = 1, r = 0. 



(5.10) 



(5.11) 



(5.12) 



Note that in the limit A ^ one recovers the boundary conditions of the linear 
wave equation. 

As an example we choose the particular solution (|5.9p obtained for A{x) = x 
(see Fig. [TT] in [Appendix A\ . Again, we read off the corresponding initial data and 
apply our numerical method to solve the IB VP. The results are displayed in Fig. [H 
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Fig. 8. The numerical and analytical residuals for the example solution to the non-linear wave 
equation with A{x) = x, see 115.91 1. Parameters: A = 1, = Zn^, Tjnin = 0.3. Left panel: Tmax = 
0.6. Right panel: Tmax = 1-2, see Fig. Illl for a plot of the solution. 
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We have solved the hyperboloidal IB VP for both a small and a large time interval 
(left and right panel in Fig. [8|). One finds a numerical accuracy of about 10 to 11 
correct digits. Hence, the solution is nearly as precise as the corresponding example 
solution of the linear wave equation with A{x) = x, cf. Fig. IH 

6. Regularized wave equation 

An important point in the preceding discussion of the wave equation is the degen- 
eracy at ^ , which is a consequence of the conformal approach. This degeneracy 
provided us with a particular first order boundary condition to be imposed at 
As an alternative to the numerical solution of this singular equation, we investigate 
in this section a regularized version. The motivation behind this study is our future 
goal of solving the Einstein equations numerically on a conformally compactified 
space-time. Note that a regular formulation of Einstein's equations can be identified 
(Friedrich's regular conformal field equations, [9]). 

For the derivation of the regularized wave equation, consider two conformally 
related n-dimensional metrics g and g with 

= ^'^gf,,y, (6.1) 
where — ft{x'^) is a conformal factor. Then, the relation 

holds, see [16]. R and R denote the Ricci scalars of the two conformal metrics, and 
□ and □ are the wave operators^. In our case we choose fl = PQ [cf . (|2.6p ] . With 
n = A, R ~ (flat Minkowski space-time), R = 6, and with the definition / i^~^f 
we obtain 

af^a{nf)^n^a~i)f. (6.3) 

Hence, equivalently to the singular equation □/ = 0, we may study the regular 
equation 

(□-l)/4a..-/>.) + ^-/ = 0. (6.4) 

This equation does not reduce to a first order boundary condition at J^"*". For the 
completeness of the algebraic system discussed in section 13.41 we therefore need to 
require the nonlinear second order wave equation as a condition there. (This second 
order boundary condition can be imposed since no "information from outside" can 
encounter the numerical domain, i.e. no ingoing characteristic crosses the boundary 
in question.) 

Again we study an explicit example by looking at the hyperboloidal IBVP (C) 
in Tab.[Tl but this time we solve the regular wave equation (j6.4p . It turns out that 

^ The wave operators are defined by □ = §'"^V^Vi/ and □ = g'^'^V^V^. 
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this formulation of the problem also permits highly accurate numerical solutions. 
However, instead of 12 correct numerical digits (as in Fig. [6p) we obtain 11 digits, 
thereby loosing one order in the numerical accuracy. 

Note that / possesses a similar accuracy as the normal derivative df /dn with 
respect to ingoing characteristics. In particular, the content of outgoing radiation 
at J^'^, described by 



is similarly precisely given by / (calculated from the regularized wave equation) 
and by df /dn (calculated from the singular wave equation). 

We conclude that, remarkably, our spectral algorithm is applicable to both reg- 
ular and degenerate wave equations. For the fully spectral scheme, the particular 
formulation of a given equation seems to play a subordinate role. 

7. Discussion 

In this paper, we constructed numerical solutions of hyperbolic equations utilizing a 
fully pseudospectral scheme. Combining this method with the concept of conformal 
infinity, we were able to obtain highly accurate solutions. 

Interestingly, by means of the method presented, hyperbolic equations can be 
handled quite similarly to elliptic equations. There is no principal difference between 
the treatment of the spatial coordinate a and that of the time coordinate t in the 
Chebyshev approximation ()3.8|1 . 

The formulation of boundary/initial values is, however, fundamentally different 
for elliptic and hyperbolic problems. In an elliptic problem, at each boundary one 
first order condition is required, i.e. a condition which contains no second order 
normal derivative with respect to the boundary in question. On the other hand, 
a hyperbolic problem requires knowledge of the unknown function and its normal 
(time-) derivative at the initial slice, i.e. two first order conditions there. The future 
boundary T = T,„ax is treated like an interior point, i.e. the hyperbolic equation 
provides a second order condition there. Note that any attempt to impose a first 
order condition at this boundary must fai|j| — not only for our pseudospectral 
algorithm but for any numerical scheme. 

The method being presented possesses a highly implicit character. The time 
evolution is not performed successively, moving from one time slice to the next, 
but the entire system is solved simultaneously instead. As a consequence, it turned 
out that for all our numerical examples it was not necessary to respect a Courant- 
Friedrichs-Lewy (CFL) condition. In fact, the spatial collocation points could be 
distributed much more densely than the time points. 




(6.5) 



'As an example, for the wave equation a particular Dirichlet type boundary value problem can be 
shown to admit infinitely many solutions. 
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We have demonstrated that for solutions which admit a rapidly converging 
Chcbyshcv expansion, the method leads to an accuracy of up to 12 or 13 correct 
digits (for a double precision code). In all examples being discussed, the numerical 
error (-Rnum) is close to the analytical error (i?an), i-e. the algorithm proves to be 
very effective. Furthermore, a geometric convergence rate is exhibited, i.e. the error 
decreases exponentially with the resolution. 

Our results encourage us to attempt to develop similar numerical techniques 
to solve the dynamical Einstein equations. As an interesting area of application 
we envision time evolution problems of perturbed axisymmetric equilibrium con- 
figurations (as e.g. rotating stars, rings or black holes with surrounding matter). 
Consequently, the method would permit a highly accurate stability analysis of such 
objects and a careful investigation of the emitted gravitational waves at J^. 

In order to apply our pseudospectral scheme to such problems, it will be nec- 
essary to cover the space-time with more than a single computational domain. (At 
least separate domains for matter and vacuum regions are required.) For the treat- 
ment of non-sphericaUy symmetric equations (i.e. higher dimensional problems), the 
computationally expensive direct matrix inversion in the Newton Raphson method 
needs to be replaced by an iterative inversion method, as already mentioned earlier. 
We believe that through the implementation of appropriate technical details our 
method becomes applicable to the solution of such physically interesting problems. 
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Appendix A. Plots of the example solutions 

Within this appendix we provide plots of aU numerical example solutions to the 
linear, inhomogeneous and non-linear wave equation that we have studied in this 
paper. 




Fig. 9. Solution to the linear wave equation with A{x) = x, cf. Figs. l4l and [5l 




Fig. 10. Example solution to the inhomogeneous wave equation, cf. Fig. [7] 




Fig. 11. Example solution to the non-linear wave equation, cf. Fig. |8] 
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Fig. 12. Numerical examples for the hyperboloidal IBVP of the linear wave equation in Tab.[T] cf. 
Fig.E] 



24 Jorg Hennig and Marcus Ansorg 



References 

[1] L. Andersson, P. T. Crusciel and H. Friedrich, On the regularity of solutions to the 

Yamabe equation and the existence of smooth hyperboloidal initial data for Einstein's 

field equations, Comm. Math. Phys. 149 (1992) 587-612. 
[2] L. Andersson and P. T. Crusciel, On hyperboloidal Cauchy data for vacuum Einstein 

equations and obstructions to smoothness of scri. Comm. Math. Phys. 161 (1994) 

533-568. 

[3] M. Ansorg, A. Kleinwachter, and R. Meinel, Highly accurate calculation of rotating 

neutron stars, Astron. Astrophys. 381 (2002) L49-L52. 
[4] M. Boyle, D. A. Brown, L. E. Kidder, A. H. Mroue, H. P. Pfeiffer, M. A. Scheel, 
G. B. Cook, and S. A. Teukolsky, High-accuracy comparison of numerical relativity 
simulations with post-Newtonian expansions, Phys. Rev. D 76 (2007) 124038. 
[5] J. Corvino, Scalar curvature deformation and a gluing construction for the Einstein 

constraint equations, Commun. Math. Phys. 214 (2000) 137-189. 
[6] S. Dain and H. Friedrich, Asymptotically flat initial data with prescribed regularity 

at infinity, Commun. Math. Phys. 222 (2001) 569-609. 
[7] J. Frauendiener, Conformal Infinity, Living Rev. Relativity 7 (2004), 1. URL (cited 

on 29 Nov 2007): http://www.livingreviews.org/lrr-2004-l 
[8] J. Frauendiener, Calculating initial data for the conformal field equations by pseu- 

dospectral methods, J. Comput. Appl. Math. 109 (1999) 457-491. 
[9] H. Friedrich, Conformal Einstein Evolution, in The conformal structure of space-time 
(Springer- Verlag, Berlin, Heidelberg, New York, 2002), pp. 1-50. 
[10] P. Grandclement and J. Novak, Spectral methods for numerical relativ- 
ity. Living Rev. Relativit y 12 (200 9), 1. URL (cited on 18 Feb 2009): 
http://www.livingreviews.org/lrr-2009-l 
[11] P. Hiibner, A scheme to numerically evolve data for the conformal Einstein equation. 

Class. Quantum Grav. 16 (1999) 2823-2843. 
[12] S. Husa, Numerical relativity with the conformal field equations, in Proceedings of the 
Spanish Relativity meeting, Madrid, 2001, Lecture Notes in Physics (Springer- Verlag, 
Heidelberg, Germany, 2002). 
[13] G. lerley, B. Spencer, and R. Worthing, Spectral methods in time for a class of 

parabolic partial differential equations, J. Comput. Phys. (1992) 88-97. 
[14] R. Penrose, The light cone at infinity, in Relativistic theories of gravitation (Pergamon 

Press, Oxford, 1964), pp. 369-373. 
[15] A. Ungor, A. Sheffer, R. B. Haber, and S.-H. Teng, Layer based solutions for con- 
strained space-time meshing, Appl. Numer. Math. 46 (2003) 425-443. 
[16] R. M. Wald, General Relativity (University of Chicago Press, Chicago, 1984). 
[17] A. Zenginoglu, A conformal approach to numerical calculations of asymptotically flat 
spacetimes, PhD thesis. Max Planck Institute for Gravitational Physics and Univer- 
sity of Potsdam, Germany, 2007, gr-qc/0711.0873. 



